Chapter 6 Does the antibiotic administration change the microbial community?

6.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point %in% c("Acclimation", "Antibiotics") ) 

label_map <- c(
  "Cold_wet" = "Cold population", 
  "Hot_dry" = "Warm-population",
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity")

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
  filter(time_point %in% c("Acclimation", "Antibiotics") ) %>% 
  ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(5.1645)  ( log )
Formula: richness ~ time_point * Population + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   434.8    446.1   -211.4    422.8       43 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.61979 -0.56608  0.06074  0.51882  2.39923 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.2193   0.4683  
Number of obs: 49, groups:  individual, 27

Fixed effects:
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                               3.7459     0.1623  23.073  < 2e-16 ***
time_pointAntibiotics                    -1.1523     0.1862  -6.190 6.02e-10 ***
PopulationHot_dry                         0.7055     0.2719   2.595  0.00946 ** 
time_pointAntibiotics:PopulationHot_dry  -0.3586     0.3037  -1.181  0.23769    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) tm_pnA PpltH_
tm_pntAntbt -0.466              
PpltnHt_dry -0.598  0.278       
tm_pntA:PH_  0.292 -0.611 -0.461

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  342.5018 353.3418 -165.2509

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.282619 7.923552

Fixed effects:  neutral ~ time_point * Population 
                                            Value Std.Error DF   t-value p-value
(Intercept)                              19.11713  2.078908 25  9.195758   0e+00
time_pointAntibiotics                   -11.46317  2.832826 20 -4.046551   6e-04
PopulationHot_dry                        24.97345  3.534826 25  7.064972   0e+00
time_pointAntibiotics:PopulationHot_dry -20.91701  4.793363 20 -4.363745   3e-04
 Correlation: 
                                        (Intr) tm_pnA PpltH_
time_pointAntibiotics                   -0.633              
PopulationHot_dry                       -0.588  0.373       
time_pointAntibiotics:PopulationHot_dry  0.374 -0.591 -0.632

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8258256 -0.5685215 -0.1867832  0.5823712  2.0651455 

Number of Observations: 49
Number of Groups: 27 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  203.6546 214.4946 -95.8273

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.7173285 1.688356

Fixed effects:  phylogenetic ~ time_point * Population 
                                            Value Std.Error DF   t-value p-value
(Intercept)                              5.365728 0.4446271 25 12.067930  0.0000
time_pointAntibiotics                   -0.761455 0.6038652 20 -1.260969  0.2218
PopulationHot_dry                        1.097291 0.7560383 25  1.451370  0.1591
time_pointAntibiotics:PopulationHot_dry -0.702255 1.0216421 20 -0.687379  0.4997
 Correlation: 
                                        (Intr) tm_pnA PpltH_
time_pointAntibiotics                   -0.631              
PopulationHot_dry                       -0.588  0.371       
time_pointAntibiotics:PopulationHot_dry  0.373 -0.591 -0.629

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.85993684 -0.51436620 -0.05130501  0.48944142  1.80209573 

Number of Observations: 49
Number of Groups: 27 

6.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.03530 0.035300 10.073    999  0.002 **
Residuals 47 0.16471 0.003505                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	1.9348458	0.1002347	6.065418	0.001
Population	1	2.1356198	0.1106358	6.694811	0.001
time_point:Population	1	0.8778553	0.0454773	2.751930	0.004
Residual	45	14.3548318	0.7436522	NA	NA
Total	48	19.3031527	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.051657 0.051657 9.9224    999  0.005 **
Residuals 47 0.244689 0.005206                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	2.0429592	0.11061666	7.254623	0.001
Population	1	2.8764490	0.15574623	10.214376	0.001
time_point:Population	1	0.8770558	0.04748846	3.114457	0.003
Residual	45	12.6723552	0.68614864	NA	NA
Total	48	18.4688192	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.67439 0.67439 55.932    999  0.001 ***
Residuals 47 0.56670 0.01206                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	1.8783556	0.23032747	16.358708	0.001
Population	1	0.7551942	0.09260332	6.577030	0.001
time_point:Population	1	0.3545686	0.04347786	3.087959	0.023
Residual	45	5.1670341	0.63359135	NA	NA
Total	48	8.1551525	1.00000000	NA	NA

dbRDA

#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta$time_pointAntibiotics" = "Antibiotics",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                  "subset_meta$time_pointAntibiotics" = "Antibiotics",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50"))+
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                  "subset_meta$time_pointAntibiotics" = "Antibiotics",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")

CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

6.3 Differences in MAGs and functional capacity

6.3.1 Cold population

6.3.1.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

subset_meta <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("Antibiotics","Acclimation"))

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% 
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 9 × 2
# Rowwise: 
  phylum               sample_count
  <chr>                       <int>
1 p__Bacillota_A                 43
2 p__Bacteroidota                15
3 p__Bacillota                    8
4 p__Pseudomonadota               7
5 p__Cyanobacteriota              3
6 p__Verrucomicrobiota            2
7 p__Bacillota_B                  1
8 p__Bacillota_C                  1
9 p__Fusobacteriota               1

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)

# A tibble: 81 × 7
# Rowwise: 
   genome                phylum               class               order                 family                 genus          species
   <chr>                 <chr>                <chr>               <chr>                 <chr>                  <chr>          <chr>  
 1 AH1_2nd_1:bin_000006  p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Kineothrix  s__    
 2 AH1_2nd_19:bin_000005 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Roseburia   s__    
 3 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia s__    
 4 AH1_2nd_2:bin_000003  p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__JAAYNV01    s__    
 5 AH1_2nd_10:bin_000007 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__                    g__            s__    
 6 AH1_2nd_15:bin_000043 p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Bacteroidaceae      g__Bacteroides s__    
 7 AH1_2nd_10:bin_000041 p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Rikenellaceae       g__Alistipes   s__    
 8 AH1_2nd_8:bin_000007  p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma      s__    
 9 AH1_2nd_14:bin_000015 p__Bacillota_A       c__Clostridia       o__Christensenellales f__UBA3700             g__            s__    
10 AH1_2nd_5:bin_000001  p__Bacillota_C       c__Negativicutes    o__Selenomonadales    f__                    g__            s__    
# ℹ 71 more rows

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 4 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus         species            
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>         <chr>              
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus    s__Proteus cibarius
2 LI1_2nd_7:bin_000001  p__Bacillota_A    c__Clostridia          o__Clostridiales      f__Clostridiaceae     g__Sarcina    s__                
3 AH1_2nd_7:bin_000055  p__Bacillota      c__Bacilli             o__Mycoplasmatales    f__Mycoplasmoidaceae  g__Ureaplasma s__                
4 AH1_2nd_13:bin_000025 p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA3700            g__           s__                

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

6.3.1.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, 
     file = "data/ancombc_antib_accl_2003.Rdata")

load("data/ancombc_antib_accl_2003.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
  filter(p_time_pointAntibiotics < 0.05) %>%
  dplyr::arrange(p_time_pointAntibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_pointAntibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

Acclimation

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

            phylum sample_count
1     Bacteroidota           15
2      Bacillota_A            5
3        Bacillota            1
4 Campylobacterota            1
5  Cyanobacteriota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  select(phylum, genus, species)

             phylum           genus species
1  Campylobacterota  Helicobacter_J        
2       Bacillota_A                        
3      Bacteroidota     Odoribacter        
4      Bacteroidota     Bacteroides        
5      Bacteroidota     Bacteroides        
6         Bacillota                        
7      Bacteroidota     Bacteroides        
8      Bacteroidota     Phocaeicola        
9      Bacteroidota     Odoribacter        
10     Bacteroidota     Phocaeicola        
11     Bacteroidota     Bacteroides        
12     Bacteroidota     Bacteroides        
13     Bacteroidota Parabacteroides        
14     Bacteroidota                        
15     Bacteroidota Parabacteroides        
16  Cyanobacteriota       Scatousia        
17      Bacillota_A                        
18      Bacillota_A    Lacrimispora        
19     Bacteroidota       Alistipes        
20      Bacillota_A                        
21      Bacillota_A                        
22     Bacteroidota       Alistipes        
23     Bacteroidota     Odoribacter        

Antibiotics

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

             phylum sample_count
1       Bacillota_A            3
2         Bacillota            1
3    Pseudomonadota            1
4 Verrucomicrobiota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  select(phylum, genus, species)

             phylum          genus species
1       Bacillota_A     MGBC116941        
2         Bacillota        CAG-345        
3       Bacillota_A      Blautia_A        
4       Bacillota_A Intestinimonas        
5 Verrucomicrobiota                       
6    Pseudomonadota                       

6.3.1.3 Differences in the functional capacity

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_filtered <- genome_counts_filt %>% 
    select(one_of(c("genome",sample_sub$sample))) 
  
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

genome_counts_filtered_row <- genome_counts_filtered %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)

6.3.1.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94332, p-value = 0.09304

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=c('#008080', '#00808050')) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

6.3.1.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95742, p-value = 0.2331

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	treatment
B06	0.68110280	0.47325490	Organic anion biosynthesis	0.20784790	Acclimation
B02	0.59963040	0.41562100	Amino acid biosynthesis	0.18400940	Acclimation
D02	0.38587770	0.20636760	Polysaccharide degradation	0.17951010	Acclimation
S03	0.27103471	0.09357224	Spore	0.17746247	Acclimation
B01	0.84215140	0.69078910	Nucleic acid biosynthesis	0.15136230	Acclimation
B07	0.44558700	0.30432910	Vitamin biosynthesis	0.14125790	Acclimation
B08	0.44596150	0.32044710	Aromatic compound biosynthesis	0.12551440	Acclimation
D09	0.25533360	0.13438350	Antibiotic degradation	0.12095010	Acclimation
B04	0.54457570	0.42921780	SCFA biosynthesis	0.11535790	Acclimation
D03	0.29210750	0.20698550	Sugar degradation	0.08512200	Acclimation
D05	0.28856110	0.22258820	Amino acid degradation	0.06597290	Acclimation
B10	0.05947806	0.03621687	Antibiotic biosynthesis	0.02326119	Acclimation

6.3.2 Warm population

6.3.2.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

sample_sub <- sample_metadata %>%
  filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Antibiotics"))

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",sample_sub$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% 
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 12 × 2
# Rowwise: 
   phylum               sample_count
   <chr>                       <int>
 1 p__Bacillota_A                 63
 2 p__Bacteroidota                32
 3 p__Bacillota                   16
 4 p__Cyanobacteriota              9
 5 p__Pseudomonadota               7
 6 p__Bacillota_B                  3
 7 p__Desulfobacterota             3
 8 p__Bacillota_C                  2
 9 p__Verrucomicrobiota            2
10 p__Actinomycetota               1
11 p__Campylobacterota             1
12 p__Elusimicrobiota              1

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)

# A tibble: 140 × 7
# Rowwise: 
   genome                phylum               class               order                 family                 genus            species
   <chr>                 <chr>                <chr>               <chr>                 <chr>                  <chr>            <chr>  
 1 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia   s__    
 2 LI1_2nd_8:bin_000013  p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia   s__    
 3 AH1_2nd_19:bin_000005 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Roseburia     s__    
 4 AH1_2nd_20:bin_000009 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Kineothrix    s__    
 5 AH1_2nd_18:bin_000011 p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Coprobacillaceae    g__Coprobacillus s__    
 6 LI1_2nd_2:bin_000001  p__Bacillota_A       c__Clostridia       o__Oscillospirales    f__Ruminococcaceae     g__              s__    
 7 AH1_2nd_6:bin_000035  p__Bacillota         c__Bacilli          o__RF39               f__UBA660              g__Faecisoma     s__    
 8 LI1_2nd_4:bin_000003  p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma        s__    
 9 AH1_2nd_1:bin_000039  p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Marinifilaceae      g__Odoribacter   s__    
10 AH1_2nd_2:bin_000018  p__Elusimicrobiota   c__Elusimicrobia    o__Elusimicrobiales   f__Elusimicrobiaceae   g__UBA1436       s__    
# ℹ 130 more rows

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 7 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus              species                       
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>              <chr>                         
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus         s__Proteus cibarius           
2 LI1_2nd_8:bin_000030  p__Actinomycetota c__Actinomycetia       o__Mycobacteriales    f__Mycobacteriaceae   g__Corynebacterium s__                           
3 LI1_2nd_8:bin_000077  p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__Parabacteroides sp003480915
4 LI1_2nd_5:bin_000032  p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA1242            g__Caccovivens     s__                           
5 AH1_2nd_18:bin_000013 p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__                           
6 LI1_2nd_5:bin_000069  p__Bacillota      c__Bacilli             o__RF39               f__UBA660             g__                s__                           
7 AH1_2nd_20:bin_000061 p__Bacillota      c__Bacilli             o__Lactobacillales    f__Enterococcaceae    g__Enterococcus    s__                           

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 20 × 3
   trait p_value p_adjust
   <chr>   <dbl>    <dbl>
 1 B0220 0.00618   0.0497
 2 B0708 0.00280   0.0497
 3 B0709 0.00662   0.0497
 4 B1012 0.00618   0.0497
 5 D0203 0.00559   0.0497
 6 D0205 0.00280   0.0497
 7 D0206 0.00280   0.0497
 8 D0207 0.00280   0.0497
 9 D0210 0.00280   0.0497
10 D0211 0.00662   0.0497
11 D0213 0.00685   0.0497
12 D0305 0.00685   0.0497
13 D0308 0.00685   0.0497
14 D0610 0.00618   0.0497
15 D0706 0.00662   0.0497
16 D0902 0.00618   0.0497
17 D0905 0.00618   0.0497
18 D0908 0.00280   0.0497
19 D0911 0.00618   0.0497
20 S0104 0.00685   0.0497

6.3.2.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, 
     file = "data/ancombc_antib_accl_warm_2003.Rdata")

load("data/ancombc_antib_accl_warm_2003.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
  filter(p_time_pointAntibiotics < 0.05) %>%
  dplyr::arrange(p_time_pointAntibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_pointAntibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + 
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs Acclimation

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  select(phylum, genus, species)

        phylum           genus                  species
1 Bacteroidota     Phocaeicola                         
2 Bacteroidota Parabacteroides Parabacteroides gordonii

Antibiotics

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

          phylum sample_count
1   Bacteroidota            3
2      Bacillota            1
3    Bacillota_A            1
4    Chlamydiota            1
5 Pseudomonadota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  select(phylum, genus, species)

          phylum       genus                  species
1   Bacteroidota Bacteroides                         
2   Bacteroidota Bacteroides Bacteroides intestinalis
3    Bacillota_A  MGBC116941                         
4      Bacillota  Ureaplasma                         
5   Bacteroidota Bacteroides                         
6 Pseudomonadota                                     
7    Chlamydiota                                     

6.3.2.3 Differences in the functional capacity

sample_sub <- sample_metadata %>%
  filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Antibiotics"))

genome_counts_filtered <- genome_counts_filt %>% 
    select(one_of(c("genome",sample_sub$sample)))

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

genome_counts_filtered_row <- genome_counts_filtered %>%
  mutate_at(vars(-genome),~./sum(.))  %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)

6.3.2.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.355 0.0228
2 Antibiotics 0.253 0.0673

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.8813, p-value = 0.03338

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 69, p-value = 0.0005759
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

6.3.2.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.348 0.0158
2 Antibiotics 0.268 0.0586

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.85785, p-value = 0.01416

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 69, p-value = 0.0005759
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	treatment
D02	0.42786820	0.2181666	Polysaccharide degradation	0.20970160	Acclimation
B02	0.62818050	0.4716670	Amino acid biosynthesis	0.15651350	Acclimation
D03	0.34025420	0.1882518	Sugar degradation	0.15200240	Acclimation
B08	0.44464280	0.3309629	Aromatic compound biosynthesis	0.11367990	Acclimation
D05	0.33295630	0.2247075	Amino acid degradation	0.10824880	Acclimation
D09	0.24723300	0.1422363	Antibiotic degradation	0.10499670	Acclimation
B06	0.70995740	0.6122709	Organic anion biosynthesis	0.09768650	Acclimation
B09	0.02951089	0.0156927	Metallophore biosynthesis	0.01381819	Acclimation